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We present an ab initio calculation of small numbers of trapped, strongly interacting fermions 
using the Green's Function Monte Carlo method (GFMC). The ground state energy, density profile 
and pairing gap are calculated for particle numbers N = 2 ~ 22 using the parameter-free "unitary" 
interaction. Trial wave functions are taken of the form of correlated pairs in a harmonic oscillator 
basis. We find that the lowest energies are obtained with a minimum explicit pair correlation 
beyond that needed to exploit the degeneracy of oscillator states. We find that energies can be 
well fitted by the expression utfEtf + Amod(iV, 2) where Etf is the Thomas-Fermi energy of a 
nonintcracting gas in the trap and A is a pairing gap. There is no evidence of a shell correction 
energy in the systematics, but the density distributions show pronounced shell effects. We find the 
value A = 0.7 ± Q.2u for the pairing gap. This is smaller than the value found for the uniform gas 
at a density corresponding to the central density of the trapped gas. 
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The physics of cold trapped atoms in quantum conden- 
sates has seen remarkable advances on the experimental 
front, particularly with the possibility to study pairing 
condensates in fermionic systems 0, 0, 0, 0, S S 13 • Many 
features of systems in the size range N ~ 10 4 — 10 6 
are now well-explored, but the small-iV limit is also of 
great interest for optical lattices. In this work we investi- 
gate the properties of small systems of trapped fermionic 
atoms using the Green's Function Monte Carlo technique 
(GFMC) that has been successful in the study of the ho- 
mogeneous gas[j|. The small systems are in some ways 
more challenging because simplifications that follow from 
translational invariance are not present. Our main goal 
here is to see how the bulk behavior evolves as a function 
of the number of atoms and to provide benchmark ab ini- 
tio results to test other theoretical methods. The Hamil- 
tonian for interacting atoms in a spherical harmonic trap 
is given by 
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where u> is the trap frequency and v(r) is the interaction 
between the atoms of opposite spin states. We will use 
units with h = 1. The interaction is chosen to approach 
the so-called unitary limit, meaning that it supports a 
two-body bound state at zero energy as well as having a 
range much shorter than any other length scales in the 
Hamiltonian. For technical reasons, we keep interaction 
range finite in the GFMC using the form 
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The effective range of the potential is ro; the short-range 
limit is roy^muj -C 1. The results are presented in the 



following sections, together with some comparison to the 
expectations based on the local density approximation 
(LDA). 



To apply the GFMC, one starts with a trial wave func- 
tion i&T that is antisymmctrized according to the fcrmion 
statistics of the particles. The GFMC gives the lowest en- 
ergy state in the space of wave functions that have the 
same nodal structure as ^t- We tried several approaches 
to parameterize ^t- For even particle number N, they 
can all be expressed in the form 
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Here <p (2 \i,j) is a pair wave function, iVf = iVj, = 
N/2, and IT/ is a Jastrow correlation factor. The an- 
tisymetrization operation A is carried out by evaluating 
determinant of an A r /2-dimensional matrix. For systems 
with an odd number of particles, we need to include an 
unpaired particle in the wave function. We define an 
orbital wave function ip{ r ) f° r the extra particle and it 
is included by adding an extra row and column to the 
determinant [8|, Eq. 14] in the antisymmetrization. We 
have mostly investigated trial wave functions where the 
pair state takes the form [nj EH , 

^)( ri ,r 2 ) = X> A £]^]>>l) ^ +'7^/2TTlx 
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Here ip n im is the oscillator state labeled by radial quan- 
tum number n and angular momentum quantum num- 
bers I, m, the oscillator shell is A = 2n+l, and A c is a shell 
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cutoff. Clebsch-Gordan coefficients (-l) l+m / s/TTl al- 
lows that the the pairs with angular momenta (I, m) and 
(I, — m ) to form zero total angular momentum state. The 
Ansatz Eq. (J4j) is analogous to the pair wave function used 
to calculation the uniform system 0] . There the particle 
orbitals were plane waves and each was paired with the 
orbital of opposite momentum. This pair state allows 
for intra-shell (n = n') as well as multi-shell (n ^ n') 
pairings. At shell closures such as N — 8, 20 the trial 
function is equivalent to the Slater determinant of har- 
monic oscillator orbitals when the cutoff is at the highest 
occupied shell and multi-shell pairings are neglected. We 
have also considered taking the pair wave function as the 
eigenstate of the two-particle Hamiltonian, requiring in 
principle an infinite cutoff in the oscillator representation. 
We call this case as 2B. 

We carry out the GFMC in the usual way described 
in Ref. @. The ground state wave function is projected 
out of the trial wave function by evolving it in imaginary 
time, and the energy is taken by the normalized matrix 
element of the Hamiltonian operator. This may be ex- 
pressed as 

|* ) = lim e- nT \^ T ) (5) 
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and E a = linv^oo {^ T He- HT \^ T ) / {^ T er nT \^ T ) . The 
integral is evaluated by the Monte Carlo method, car- 
rying out the exponentiation by the expansion e~ HT w 

( e -VAT/2 e -TAT e -VAr/2)M and usmg path sampling . 

Our target accuracy is 1% on the energies. This is 
achieved by taking numerical parameters At = 0.04w 
and 15, 000 < M < 30, 000. In practice, the convergence 
to the ground state is reached in the first few thousands 
of steps. The Monte Carlo sample points that leave the 
region where \^t) > are discarded. This nodal con- 
straint avoids the signal decay known as 'fermion sign 
problem'. The energies depend on the range parameter 
ro only in the combination rg y/mu; which we set to 0.1. 
We believe this is small enough to give energies that ap- 
proach the contact limit to within 1%. Smaller values of 
range parameter are possible but increase the statistical 
fluctuations of the Monte Carlo integration. 

The cases N = 2, 3 are special in that analytic solu- 
tions are known. For the N = 2 system, the Jastrow 
correlation factor can be defined to give the exact wave 
function and energy E — 2u). The TV = 3 system gives 
the first real test of the theory. The exact energy is 
E = 4.27.. .a;, given by the solution of a transcendental 
equation in one variable 12). Using Eq. [4] with a single 
term (A c = 0) we find an energy of 4.28 ± 0.04w in close 
agreement with the exact value. In contrast, taking the 
pair wave function as the two-particle eigenstate gave a 
significant difference, 4.41 ± 0.02w. 



We now turn to the larger systems and determine the 
parameters in Eq. |4j As mentioned earlier, at the 
shell closures taking the cutoff at the highest occupied 
shell gives the harmonic oscillator Slater determinant 
(HOSD). We also tried Slater determinants for mid-shell 
systems, but typically they break rotational symmetry 
and give a significantly higher energy. Use of Eq. [4] guar- 
antees that will be rotational invariant (for even N). 
One parameterization we explored was to use the results 
of Ref. [H to guide the choice of A c and the «a- This 
gives a rather open trial wave function, having signifi- 
cant particle excitation out of the lower shells. Another 
choice is to take each qa proportional to the shell oc- 
cupancy of the HOSD, which we call II. For mid-shell 
systems, this also produces significant excitations out of 
the nominally closed shells. Both these schemes gave 
poorer energies than we could achieve by taking ampli- 
tudes «a that maximize the occupancy of the filled shells 
of the HOSD, and have no occupancy in the nominally 
empty shells. Values of a a that approach this compact 
limit (CL) are given in Table [J As seen in this table, a a 
parameters are not very sensitive in the ranges limited 
by the shell closures and are kept constant. 

For odd N, we take the pair wave function the same 
as in the neighboring even system, and the orbital of 
the odd particle as an oscillator state of the filling shell. 
Thus, in the range N = 3 — 7 the orbital is a p-shell 
orbital with (n,l) — (0, 1). Starting at A = 2 there is a 
choice of orbitals, eg. (n,l) = (0,2) or (1,0) for N = 9. 
We found for the N = 9, 11, 13, 15 and 19 systems, the 
energies are degenerate within the statistical errors and 
it was not possible to determine the density preference of 
the excitation. For these cases, we simply took the odd 
orbital to be one with the highest value of I. 

The calculated energies are summarized in Table |TT] 
for the paired wave function in the compact limit. The 
statistical errors of the GFMC are given in parenthesis. 
One sees that an accuracy of 1% is achieved with the 
numerical procedure we described earlier. In Fig. [T] we 
show a plot of the energies including the results from the 
HOSD trial wave function. As expected, the energies are 
the same at the shell closures but the HOSD gives higher 
energies in mid-shell. As an example of the sensitivity 
of the energy to the detailed assumption about the pair 
state, the results for N = 12 are: E C l = 21.5(3), E 2 b = 
22.3(2) and E n = 22.4(3). One sees that the energies 
are actually rather close. However, Eql is consistently 
2—4% below other pairing node assumptions in the range 
of N considered. In case of the simple HOSD, Ehosd = 
23.0(1) which is 7% above the CL energy. This gives 
some confidence that the assumed nodal structure of $>x 
is adequate for our purposes. We will comment on the 
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TABLE I: Parametrization of qa (CL). 



QO Ql 


02 


«3 




1.0 








N = 2,3 


1.0 0.1 








4 < N < 7 


1.0 1.0 








N = 8,9 


1.0 0.5 


0.01 





10 < N < 19 


1.0 1.0 


1.0 





N = 20,21 


1.0 1.0 


0.5 


0.01 


N > 22 



TABLE II: GFMC energies for the unitary trapped fermion 
gas with the CL pair functions. Also shown are the energies 
of noninteracting gas (HOSD). The unit of energy is u. 



N 


HOSD GFMC 


N 


HOSD GFMC 


2 


3 


2.01(2) 


13 


35.5 


25.2(3) 


3 


5.5 


4.28(4) 


14 


39 


26.6(4) 


4 


8 


5.1(1) 


15 


42.5 


30.0(1) 


5 


10.5 


7.6(1) 


16 


46 


31.9(3) 


6 


13 


8.7(1) 


17 


49.5 


35.4(2) 


7 


15.5 


11.3(1) 


18 


53 


37.4(3) 


8 


18 


12.6(1) 


19 


56.5 


41.1(3) 


9 


21.5 


15.6(2) 


20 


60 


43.2(4) 


10 


25 


17.2(2) 


21 


64.5 


46.9(2) 


11 


28.5 


19.9(2) 


22 


69 


49.3(1) 


12 


32 


21.5(3) 









nodal structure again later. 

These results show that the pairing is less important 
in the trial wave function for the finite systems than it is 
in the uniform gas. While in the homogeneous gas a BCS 
treatment of the trial function lowers the energy by more 
than 20% {E SF = 0.42£ FG and E norrnal = 0.56£ FG ), 
the difference in energy between the both phases of the 
trapped gas does not exceed 7% at the most in the open 
shell configuration N = 12. At the shell closures, the 
BCS treatment does not offer any improvement 

There is a virial theorem for unitary trapped gases 
given by[I3 E„ = 2N(U), where (U) = \muj 2 (r 2 ). 
The theory can thus be tested by independently cal- 
culating the expectation value of the trapping poten- 
tial. Expectation values of operators are usually esti- 
mated in the GFMC by the expression (* |W|*o) ~ 
2(U)gfmc — (U)var.- Using this estimate of (U), we find 
energies somewhat lower than obtained by direct GFMC 
calculation (see Fig. [1]). This could be due to the errors 
associated with the extrapolation formula for expectation 
values. 



We now examine how well the energies fit the asymp- 
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FIG. 1: (color online) Energy systematics of the trapped uni- 
tary Fermi Gas. Circles, GFMC calculated with the HOSD 
trial function; triangles, GFMC with the CL trial function; 
crosses, virial formula. The dotted line is the energy of the 
HOSD for free particles. The dashed line is the fit (Eq. [6} to 
the CL calculated energies. The unit of energies is cu. 



totic theory for large nonuniform systems. The first term 
in the theory is the Thomas-Fermi (TF) approximation 
[l4T | ; the TF approximation to the trapped unitary Fermi 
gas is Etf(0 = £ 1/2 w(3A0 4/3 /4, where £ is the univer- 
sal constant relating the energy of the uniform gas to 
that of the free Fermi gas. Adding a second term in the 
expansion gives a better description of the energy of the 
harmonic oscillator energy of the trapped gas in the large 
N limit [lij ]. We therefore will include that in the fit to 
the energies, using the form 



E' TF (0 



1/2 



U> 



/ (3AQ 4 / 3 (3AQ 2 / 3 
I 4 + 8 



(6) 



As it may be seen from Fig. [T] there is also a significant 
odd-even variation in the energies. We shall include this 
effect as well by fitting to the function 



E = E' TF {£) + Amod(AT,2) 



(7) 



The result of the fit is shown by the dashed line in Fig. 
[TJ The fit value of £ is £ = 0.50. This is somewhat higher 
than the bulk value 0.42 — 0.44. This suggests that the 
convergence to the bulk is rather slow. One might expect 
to see shell effects in the energies once the smooth trends 
have been taken out. The HO energies, for example, os- 
cillate around Etf(0) with (negative) peaks at the shell 
closures. The effect is visible in the abrupt change of 
slop of the HO curve in Fig. [T] However, in our fit to the 
calculated energies, we do not see a visible shell closure 
effect. In the fit we find for the parameter A the value 
A = 0.6w, in accordance with the average of the odd-even 
staggering of the energy A = 0.7(2)w. If the pairing gap 
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FIG. 2: (color online) Radial densities for N = 2, 8, 14 and 
20. For N = 20, free particle density distribution is also 
shown (dotted line). The unit on the radial axis is } . 



were controlled by the density at the center, it would be 
much larger; of the order of shell spacing or higher. On 
the other hand, if the odd particle is most sensitive to the 
surface region, the pair effect could be smaller. A more 
systematic approach to LDA exists where the superfluid 



correlation is introduced ah initio 



We now turn to the density distributions calculated 
with the GFMC. We determine the density by binning 
the values obtained by the Monte Carlo sampling. With 
~ 15000 paths and 1000 samples per path there is ade- 
quate statistics to get details of the density distribution 
well beyond the mean square radii. Fig. [2] shows the cal- 
culated densities for N = 2, 8, 14, and 20. We notice that 
the central densities show a pronounced dip for N = 8 
and a peak for N = 20. These are characteristic of shell 
structure, depending on whether the highest occupied 
shell has s-wave orbitals or not. Fig. [5] also shows the 
HO density for N — 20. The density of the interacting 
system is more compact, as required by its lower energy 
and the virial theorem. The central peak has roughly the 
same relative shape in the two cases. Thus the basic HO 
pattern is maintained, even though the system shrinks in 
size. 



Let us return again to the problem of the trial func- 
tion and its fixed nodal structure. It is interesting to 
ask how different are the nodal positions for the different 
fr's. We can characterize the trial wave functions by 
their relative overlaps of the sign domains. If we define 



the nodal overlap between and 4'r2 is given as 

Nodal overlap = max[100% — x,x\. From this defini- 
tion, nodal overlap ranges 50% ~ 100%. Because of 
the strong suppression of the superfluidity at the shell 
closures the nodal overlap between the normal and the 
superfluid node wave functions seems to be the largest 
~ 100% at the shell closures and has low values ~ 85% 
at N = 5 and - 67% at N = 14. 



We believe our computed energy systematics is reli- 
able enough to arrive at the following conclusions. 1) 
The energies are significantly higher than given by the 
TF model with bulk £ = 42 - 44. 2) Stabilization of 
closed shell systems with respect to open shell ones is 
much weaker than in the free gas. However, the den- 
sity distribution has pronounced fluctuations similar to 
those of the pure harmonic oscillator density. 3) There is 
a substantial pairing visible in the odd-even binding en- 
ergy differences, but the magnitude is less than the bulk 
pairing parameter associated with a uniform system of 
density equal to the central value in the finite system. 
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